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We simulate crystallization and melting with local Monte Carlo (LMC), event-chain Monte Carlo (ECMC), 
and with event-driven molecular dynamics (EDMD) in systems with up to one million three-dimensional hard 
spheres. We illustrate that our implementations of the three algorithms rigorously coincide in their equilibrium 
properties. We then study nucleation in the NVE ensemble from the fee crystal into the homogeneous liquid 
phase and from the liquid into the homogeneous crystal. ECMC and EDMD both approach equilibrium orders 
of magnitude faster than LMC. ECMC is also notably faster than EDMD, especially for the equilibration 
into a crystal from a disordered initial condition at high density. ECMC can be trivially implemented for 
hard-sphere and for soft-sphere potentials, and we suggest possible applications of this algorithm for studying 
jamming and the physics of glasses, as well as disordered systems. 
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I. INTRODUCTION 

Crystallization and melting have long been central sub¬ 
jects in statistical physics. These processes connect mi¬ 
croscopic nucleation with the macroscopic phenomena of 
domain growth and of phase transitions. A number of 
numerical methods and simulation techniques have been 
brought to bear on these subjects, following the pioneer¬ 
ing computer simulations of hard-sphere systems by both 
Monte Carlo 1 and by molecular dynamicJ®^. 

The hard-sphere system is trivial to describe. Never¬ 
theless, equilibration in this simplest of all particle sys¬ 
tems is a slow process, because of the large activation 
free energy for crystallization. Timescales are also espe¬ 
cially large in the fluid-solid coexistence regime, because 
of the surface tension between coexisting phases. Spe¬ 
cialized algorithms for equilibration have been developed 
to overcome these problems, and the melting and crystal¬ 
lization time scales provide useful benchmarks for their 
comparison. 

A rejection-free hard-sphere “event-chain” Monte 
Carlo algorithm (ECMC) 9 has recently allowed to speed 
up equilibration for two-dimensional hard disks by 
roughly two orders of magnitude compared to the event- 
driven mole cular dynamic d 5 * 10 * (EDMD) and to local 
Monte Carlo 11 ^ (LMC). In ECMC, a randomly sam¬ 
pled starting sphere moves along a straight line until the 
latter collides with another sphere, which then moves in 
the same direction until it collides itself with yet another 
sphere. This continues until the spheres’ total displace¬ 
ment equals a certain fixed length L c . ECMC breaks 
detailed balance (moves are in the +x, -\-y, and +z direc¬ 
tions only) yet satisfies global balance and ergodicity 13 . 
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It rigorously samples the equilibrium Boltzmann distri¬ 
bution. Considerable speedup was also demonstrated for 
the extension of ECMC to continuous potentials 13 — . 

In this paper, we assess the speed of ECMC, EDMD, 
and LMC not by computing autocorrelation functions in 
equilibrium, but rather by the time scales associated with 
melting and crystallization in systems of many spheres 
at high density. We expect our observations to extend to 
subjects as dense packing, nucleation and jamming. We 
focus on the melting from the metastable solid branch 
to the stable liquid (that is, slightly below the liquid- 
solid coexistence interval) and also the nucleation pro¬ 
cesses from the metastable liquid branch towards the 
stable solid slightly above coexistence. The paper is or¬ 
ganized as follows: Our model system and our methods 
are described in Section [IIJ together with the observables 
on which we focus: The non-dimensional pressure and 
the local and global orientational order parameters. Re¬ 
sults are summarized in Section m We reproduce the 
phase diagram by ECMC and quantify the efficiencies of 
our three methods. We discuss the relative efficiency of 
ECMC and EDMD for the crystallization process. Con¬ 
cluding remarks are described in Section El 


II. MODEL, ALGORITHMS, AND OBSERVABLES 

We consider N monodisperse hard spheres of radius a 
in a cubic box of sides L with periodic boundary condi¬ 
tions (PBC). The density (packing fraction) v is given 
by v = 4/3A/brcr 3 /L 3 . We concentrate on melting and 
crystallization from the unstable to the stable phase, i.e., 
from the unstable crystalline branch to the liquid and 
from the unstable liquid branch to the crystal. For our 
melting runs, at density v = 0.490 below the coexistence 
interval of liquid and solid phases, we prepare the ini¬ 
tial configurations as perfect fee crystals, corresponding 
to the stable phase at high densities (the free energy dif- 








Hard-sphere melting and crystallization with event-chain Monte Carlo 


2 




FIG. 1. Pressure P* (left) and local order parameter q§ (right) in the dense liquid at v — 0.490, obtained from ECMC and 
EDMD for N = 131,072 hard spheres. In ECMC, the pressure is computed using the excess displacement method of Eq. <§• 
In the right panel, the local order parameter of LMC is also shown. 


ference between fee and hep c rystal have been discussed 
actively since the works of Ref) 15 l 16 l). The fee initial con¬ 
ditions are compatible with the cubic simulation box. 

For the crystallization runs at density v = 0.548 above 
the coexistence interval, we start from a fluid initial con¬ 
figuration at a liquid-phase density v = 0.490. In order 
to reach the higher target density, we repeatedly increase 
a slightly and remove all created overlaps by sliding over¬ 
lapping pairs of spheres for a half length of overlap along 
their common symmetry axis. This is done until all pair 
overlaps have disappeared. 

In EDMD, hard spheres evolve in continuous physical 
time through collisions, and the dynamics solves New¬ 
ton’s classical equations of motion. We use an efficient 
sequential implementation 10 . LMC and ECMC are im¬ 
plemented very simply. For the former, the optimal dis¬ 
placement of spheres is determined by short LMC runs 
from the initial conditions so that the acceptance ratio is 
1/2. For the latter, a single parameter, the chain length 
P c , must be optimized for each density. A single event 
can be implemented very quickly in ECMC, as the mo¬ 
tion is always in +x or +y, which decreases the CPU 
time per event. We consider systems with N = 2, 048, 
131,072 and 1,048,576 spheres. Our calculations are 
mainly done on the Intel Xeon CPU E5-2680 2.80GHz, 
where we reach ^ 3.15 x 10 9 events/h of CPU time for 
ECMC and ~ 4.62 x 10 8 collisions/h for EDMD. For the 
LMC algorithm, ^ 6.5 x 10 9 trials/h are reached. All our 
comparisons of algorithms are in terms of CPU time. To 
be as fair as possible, the three algorithms were imple¬ 
mented following unified design principles. Furthermore, 
we used the same computer, the same language (Intel 
FORTRAN), and the same optimal option of compiler. 
An event-count would have produced similar results. 

We track the time-evolution of the order¬ 


ing/disordering of the system from the pressure and 
the local and global orientational order parameters. In 
EDMD, the nondimensional virial pressure is computed 
from the collision rate via the virial theorem: 


P* = /3P(2a) 3 



/ 3m 1 
3T N 


N bij 

collisions 


(i) 


where T is the total simulation time, and f3 = 1/m (v%/) 
is the inverse kinetic temperature (mass m and mean- 
square x-component of velocity of spheres). The colli¬ 
sion force bij = r^ • v^- is defined between the relative 
positions and the relative velocities of the collision part¬ 
ners 11 17 1 In ECMC, the pressure P* can be evaluated 
from the mean excess chain displacement 13 : 


p* 


6z/ 
7r 


^final ^initial 


L c 


^ chains 


( 2 ) 


where Xfi na i and ^initial are final and initial positions of 
each chain, respectively, taking into account the PBC. 
This convenient formula replaces the tedious extrapola¬ 
tion of the pair correlation function at contact #2 (r = 2cr) 
that was used previously and that must still be used for 
LMC. The comparison of the evolving pressure of EDMD 
(using Eq. 0 ) and ECMC (using Eq. © in the liquid 
state at v = 0.490 is shown in the left of Fig.[]] The pres¬ 
sures fluctuate around P* = 11.3893, and agree within 
very tight error bars. 

Besides the pressure, we quantify the speed of melting 
and of crystallization via the time-dependent local q§ and 
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global Qq order parameters 18 : 
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where n(i) is the number of nearest neighbors for each 
sphere i and l6,m( r r/) are the spherical harmonics with 
icosahedral symmetry for the distance vector r^- be¬ 
tween spheres i and j. We detect nearest neighbors 
by the SANN algorithm of Meel et alP^ rather than by 
the Voronoi construction. On the right of Fig. [l] we 
again demonstrate that the equilibrium values for g 6 (^ 
0.38678) agree within tight error bars for LMC, ECMC 
and EDMD. Similar agreement was reached for the global 
order parameter Qq. We note that perfect fee configura¬ 
tions have an orientational order Q 6 = q 6 ~ 0.575, while 
for liquid configurations (including our disordered initial 
configurations) Qq approaches zero. In the dense liquid, 
the local order q 6 is non-zero because of the build-up of 
transient local crystal structures 20 . 


III. RESULTS 

A. Hard-sphere phase diagram 


20 


- 1 - r~ 

FCC Initial Condition 
A 1X1=1,048,576 
• 1X1=131,072 
■ 1X1=2,048 
Fluid Initial Condition 
A 1X1=1,048,576 
O 1X1=131,072 
□ 1X1=2,048 


Crystallization Simulation 

y- 


Carnahan&Starling (1969) 



15 






a** 

••• i 


10 



A 

Melting Simulation 


/O ^q' nD «^r 

. 

__ 


■ Alder&Hoover&Young(1968) 
Fernandez et al. (2012) 


0.48 


0.5 


0.52 


0.54 


FIG. 2. Hard-sphere equation of state obtained by ECMC 
after ~ (D(10 12 ) collisions. The snapshots in the coexistence 
interval at N = 1,048, 576 represent the local orientational or¬ 
der qQ (i) for each sphere i (liquid-like local order is represented 
in green - solid-like local order in blue). The equilibrium co¬ 
existence pressure for the infinite system is also shown. The 
densities v — 0.490 and v — 0.548, on which we concentrate 
in Sections |III B| and |III C[ are indicated. 


The fluid-solid coexistence in the NVE ensemble 
(which, for hard spheres, corresponds to the common 
NVT ensemble) for densities v in the interval 0.494 < 
v < 0.5 45 has been investigated for more than 50 
years 2 4 7 . Recently, various theoretical equations of 
states were compared with the results of a large-scale 
EDMD simulation on this system,^ with N ~ 10 6 . The 
metastable fluid branch in the fluid-solid coexistence win¬ 
dow was found stable against freezing on EDMD time 
scales up to ^ 0(1O 9 ) collisions. To speed up the simula¬ 
tion, the replica exchange MC method was adapted to the 
hard-sphere case. To keep the acceptance rate at reason¬ 
able values, many replicas at finely spaced densities had 
to be used, and this approach proved restricted to quite 
smal l sy stem sizes (TV = 32 and N = 108) P 2 Fernandez 
et al! 23 explored the coexistence of hard-sphere systems 
in equilibrium by tethered MC for relatively small system 
sizes ~ 0(1O 3 ). In this method, the approach to equi¬ 
librium is accelerated by a biased field of two order pa¬ 
rameters. The equilibrium pressure P* = 11.5727(10) is 
obtained through extrapolation towards the infinite-size 
limit. We note that in three-dimensional hard spheres, 
the direct simulation remains difficult in the coexistence 
region, even for ECMC, whereas for the analogous two- 
dimensional hard disks, the equilibration of the coexist¬ 
ing hexatic and liquid phases by ECMC proved possible 
at all densities, for up to one million disks. 24 


Fig-! shows the hard-sphere equation of state from 



CPU time [h] 

FIG. 3. Melting at v — 0.490 (N — 131,072) from an fee 
initial configuration into the stable liquid, tracked by the time 
evolution of the global Qq order parameter in LMC, EDMD, 
and ECMC with optimal chain length (L c /2a = 5.87). Data 
averaged over 5 samples. Note that ECMC and EDMD are 
orders of magnitude faster than LMC. 
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FIG. 4. Crystallization at v — 0.548 from a random initial configuration, tracked by the evolution of the pressure P* in 
EDMD and ECMC with different chain lengths L c . Data averaged over 100 samples for N = 131,072 (left) and 5 samples for 
N = 1,048, 576 (right). The inset in the left panel illustrates the influence of the parameter L c on the performance of ECMC. 


ECMC (final pressures after long runs with (2 ^ 3 x 10 12 ) 
collisions). The pressure is averaged over 10 10 colli¬ 
sions at the end of the simulation. The stable and 
unstable fluid pressures well agree with Hoover-Re^l 
and Carnahan-Starling extrapolation.^ Inside the coex¬ 
istence interval, the final pressure depends on the initial 
configuration, as the metastable fee solid or fluid initial 
configurations are preserved on the time-scale of the sim¬ 
ulation. In the NVE ensemble, the presence of inter¬ 
faces of different topologies makes that the equation of 
state is non-monotonous, and the liquid-solid coexistence 
pressure curve is not flat in a finite system. As one in¬ 
creases the density from the liquid phase, the spherical 
or cylindrical droplets that can be seen in Fig. [2] gen¬ 
erate an excess Laplace pressure. This is a nalog ous to 
what was found in two-dimensional hard disk J 24 * 2 ^ , which 
show droplets and stripes or in fluid-gas mixtures of the 
three-dimensional Lennard-Jones system 28 -^, where spher¬ 
ical and cylindrical droplets as well as two-dimensional 
stripes are found. 

Specifically, for v < 0.498, simulations from arbitrary 
initial conditions converge to the same pressure since 
the system is completely liquid and nucleation barriers 
are low. In the region of v = 0.500 ^ 0.514, simu¬ 
lations from fee initial configurations successfully cre¬ 
ate interfaces with different topologies, whereas simu¬ 
lations from fluid initial conditions remain on the fluid 
branch. The pressure at v — 0.498 ~ 0.514 decreases as 
P* = 12.2 ~ 11.5, and agrees with the expected coex¬ 
istence pressure.^ 3 The phase coexistence at equilibrium 
can be seen clearly by the spatial distributions of the 
local #6 order parameter, where the fee crystal reduces 
to a droplet [y = 0.500 ~ 0.502) when started from an 
fee crystal. For larger densities [y = 0.504 ^ 0.512) the 
remaining fee phase has the form of a cylinder that recon¬ 


nects through the PBC, surrounded by the liquid dom¬ 
inant phase created through melting (see the insets of 
Fig.|2|. In the density interval v = 0.514 ~ 0.530, the fee 
crystal and the liquid remain metastable on the available 
scales of simulation time. For v = 0.532 ~ 0.543, simu¬ 
lations from fluid initial conditions nucleate fee droplets. 
At v > 0.543, simulations from fluid initial condition 
drop down near the solid branch, however, relaxation is 
in progress at the value around slightly higher pressure 
than the solid branch, except for the case of N = 2, 048. 


B. Melting from an fee initial configuration at a fluid 
density 

We now study the speed of melting into the equilibrium 
liquid phase at v = 0.490 from an fee crystalline initial 
configuration for N = 131,072. From the global order 
parameter Q 6 , as shown in Fig. [3j the initial fee con¬ 
figuration (Qq ~ 0.575) rapidly becomes unstable for the 
three algorithms and approaches the liquid branch, where 
the global orientational order approaches zero. Melting 
is much faster for ECMC and EDMD than for LMC. 
ECMC is found to be somewhat faster than EDMD. 


C. Crystallization from a fluid initial condition at a solid 
density 

In Fig.[4j we show the evolution toward crystallization 
of EDMD and of ECMC with different chain lengths L c . 
Results are averaged over 100 samples for N = 131,072 
(left of Fig. [4 ) and over 5 samples for N = 1,048,576 
(right of Fig. 4b. Results for three trial runs in which 
we changed the chain length L c are also shown Fig. [4] 
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FIG. 5. Crystallization at v — 0.548 from a fluid initial configuration with N — 131,072, tracked by the evolution of local qs 
and global Qq order parameters by LMC, EDMD and ECMC with different chain lengths L c . Data averaged over 100 samples. 


The efficiency of ECMC naturally depends on L c . For 
both methods, the pressure remains somewhat above the 
configurational equilibrium pressure P* ^ 11.934. The 
relative advantage of ECMC with optimal chain length 
is evident. 

In Fig. [5] and Fig. [6j we show the evolution of the lo¬ 
cal q§ and global Qq order parameters in crystallization 
runs of EDMD and ECMC with three chain lengths. The 
number of samples is again 100 for N = 131, 072 and 5 for 
N = 1,048,576. In the early stage of relaxation, q 6 and 
Qq are increasing functions of CPU time. In the perfect 
fee configuration, the local and global order parameters 
agree to q§ = Qq = 0.57452. Due to thermal fluctua¬ 
tions, actual numerical simulation at fee equilibrium es¬ 
timate the order parameters to (qQ,Qo) = (0.505,0.483) 
in (TV, z/) = (1,048,576,0.548). Although the crystalliza¬ 
tion still proceeds, our final averaged q$ and Qq reach 
around 0.488 and 0.432, respectively. The inconsistency 
between the fee crystal structure and the (finite) simu¬ 
lation box may well be responsible for the reduction in 
order parameter. 

Order parameters decay towards higher order in time, 
ECMC with optimal chain length L c /(2a) = 23.46 needs 
48.6 CPU hours to reach at qQ = 0.4863 in TV = 131, 072, 
however, EDMD needs 165 CPU hours. Both smaller 
and larger L c results in the inefficiency of relaxation. Af¬ 
ter 165 hours, one quarter of ECMC simulations had al¬ 
most reached the equilibrium state (i.e, qQ > 0.99 x qQ eq ) 
whereas for EDMD, only 10 % had reached such values. 
The results for LMC are also shown in Fig. [5j and they 
show that it is much slower than ECMC and EDMD. 
For example, our LMC simulation on average reaches 
qQ = 48.25 after 0(1O 7 ) LMC sweeps corresponding to 
152 CPU hours. At this qQ = 48.25, ECMC and EDMD 
needs only about 11 and 28 CPU hours, respectively. The 
decay of Qq has the same tendency as that of qQ , how¬ 


ever, its values are lower than that of qQ. Since the whole 
system is slower to order than to establish local order, 
the global orientation also grows more slowly than the 
local orientation. In the larger case N = 1,048,576, 
ECMC with optimized chain length and EDMD need 
416 and 1,000 CPU hours to reach qQ = 0.4897, respec¬ 
tively. Note that if chain length is not optimized, the 
performance of ECMC changes drastically and becomes 
comparable to or slower than that of EDMD. As for qQ 
and Qq , ECMC with optimal length is faster than that 
of EDMD for a certain factor depending on the target 
point, which will be discussed in Section [ill D 

D. Relative Speed of ECMC and EDMD 

To further quantify the equilibration speed of ECMC 
and EDMD, rather than the observables as a function 
of time 0(£), we consider the elapsed CPU time T(0) 
from the beginning of the simulation t = 0 at which the 
observable O is reached. This allows us to define the 
relative efficiency R s \ 

R s (0) = 7edmd(0)/7ecmc(C ) ) 5 (5) 

and analogously for any pair of algorithms, where O is 
an observable, in our case the pressure P*, or the local 
and global order parameters qQ and Qq. 

For the melting case (not shown) at v = 0.490, R s 
for observables takes around 3.2 (t = 0) to 1 at the 
end of simulation. EDMD is slightly slower than that 
of ECMC. This relative speed remains rather unchanged 
during the melting process. The situation changes dras¬ 
tically for the crystallization process. Fig. [7[ shows R s 
as a function of normalized observables as O = (O — 
C^initiai)/(^equii. initial) v — 0.548 in the crystal¬ 
lization process, at which each data can be estimated by 
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FIG. 6. Crystallization at v — 0.548 from a random initial configuration with N — 1, 048, 576 tracked by the evolution of local 
qq and global Qq order parameters in LMC, EDMD and ECMC with different chain lengths L c . Data averaged over 5 samples. 
Note that LMC is much slower than the other two methods. 


Fig. [4]- Fig. i ^initial and Oequii. are the values at t = 0 
and at the equilibrium, respectively. Those are obtained 
by independent runs at the perfect fee crystal which are 
#6 = 0.505, Qq = 0.483, and P* = 11.934. In all cases, 
R s (0) is increasing function and growing drastically near 
the crystal (i.e., O > 0.8) as a hockey stick curve. The 
relative efficiency R s depends rather weakly on system 
size. We did not compute variations precisely, as the 
running times T were only averaged over 5 samples for 
N = 1,048,576. The different algorithmic complexity of 
our methods might marginally contribute to the size de¬ 
pendence of R s \ Our EDMD algorithm is implemented in 
0(N log N) per N collisions and ECMC as O(N) per N 
collisions. At q 6 = 0.47, R s {q6) takes 1.33 (N = 131, 072) 
and 1.49 (N = 1,048,576), respectively. 


IV. CONCLUSION 

In this paper, we compared hard-sphere Monte Carlo 
and Molecular Dynamics algorithms, that coincide in 
their equilibrium properties. In large systems with up 
to one million spheres, we recovered the known phase 
diagram and especially the coexistence region. We quan¬ 
tified the approach towards equilibrium, namely towards 
the fee crystal from the liquid-like initial configuration 
at packing v = 0.548 or the stable liquid from an fee 
initial configuration at packing v = 0.490. We clearly 
showed that the EDMD and ECMC are orders of magni¬ 
tude faster than the LMC algorithm for both the melt¬ 
ing and the crystallization. ECMC needs optimization 
for chain length L c , and we generally find that the indi¬ 
vidual chains should wrap a few times around the sim¬ 
ulation box. The effect of the chain length is rather 
drastic, and L c must be optimized carefully. The op¬ 


timal chain length for the crystallization process is es¬ 
timated around L c /(2a ) ~ 25 (N = 131,072) and 50 
(N = 1,048,576). With a fixed L c /( 2cr), the actual 
chain length (xfi na i — ^initial) /(2(7) can be obtained by 
trial and error before the production runs. It may also 
be estimated by the ECMC pressure formula, Eq. ©). 
as 


(^final ^initial) _ R (-^ / c/(2ct))7T . . 

2cr “ ^ j 

which is evolving during simulation according to pressure 
relaxation. In case of the optimal chain length L c /2a = 
23.46 ^ L^/2, the chain winds around 6 times around 
the periodic box. While doing so, very few spheres get 
hit more than once. 

We conclude that ECMC with well-chosen chain 
lengths is far superi or to LMC, although it can be imple¬ 
mented just as easil}P2En] Even with respect to molecular 
dynamics, it performs very well. The clearest advantage 
of ECMC over EDMD shows up in the crystallization, 
that is, in the buildup of long-range correlations. We 
expect the ECMC algorithm and its extension to con¬ 
tinuous potentials to be helpful to inv estiga te jamming 
and to estimate accurate nucleation rat^ 31 32 and to ana¬ 
lyze the full scenario of nucleation and precursor crystal¬ 
lization 33 . Of particular interest might be that ECMC 
remains event-driven even for continuous potentials and 
very simple to implement. Molecular dynamics, on the 
other hand, must be implemented as a time-driven al¬ 
gorithm for continuous potentials. The discretization of 
the equations of motion then makes molecular dynamics 
rather awkward to implement. 
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FIG. 7. Relative speed R s (see Eq. of LMC, ECMC and EDMD as a function of normalized observables O — (O — 
^initial)/(C>equii. — ^initial), (left) O = P* and (right) O — Qq for N — 1,048,576. As also shown in Fig. [6] both ECMC and 
EDMD are orders of magnitude faster than LMC. ECMC shows considerable advantage over EDMD in the later times of the 
evolution, when large-scale structures are built up. 
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